/****************************************************************
Compile: gcc -O3 -o design design.c
Usage: ./design mode mode1 mode2 pnflag N L inputfile
       ./design 0 0 0 1 11 14 input-2

mode          : 0 normal ; 1 cycled ; 2 triggled cycled
mode1         : 0 disallow self-activation ; 1 allow self_activation
mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II
pnflag        : 1 print the solutions
N             : number of the nodes
L             : number of the steps        

This program is used to calculate the designability for a given trajectory
under the majority rules.
****************************************************************/

#include <stdio.h>
#include <stdlib.h>

int minlink, minlinks;

void print_help ()
{
	printf ("Usage         : ./design mode mode1 mode2 pnflag N L inputfile\n");
	printf ("                (./design 0 0 0 0 11 14 input-2)\n");
	printf ("mode          : 0 normal ; 1 cycled ; 2 triggled cycled\n");
	printf ("mode1         : 0 disallow self-activation ; 1 allow self_activation\n");
	printf ("mode2         : 0 no Restrictions ; 1 Restriction I ; 2 Restriction II; 3 Restriction I and II\n");
	printf ("pnflag        : 1 print the solutions\n");
	printf ("N             : number of the nodes\n");
	printf ("L             : number of the steps\n\n");
}

void calc_links (int N, int solution[N])
{
	int i, links;

	links = 0;
	for (i = 0; i < N; i++)
	{
		if (solution[i] != 0)
			links ++;
	}
	if (links < minlink)
		minlink = links;
}

void print_solution (int N, int solution[N], int i)
{
	int j;

	printf ("Node %2d :", i);
	for (j = 0; j < N; j++)
	{
		if (solution[j] > 0)
			printf (" g");
		else if (solution[j] < 0)
			printf (" r");
		else
			printf (" n");
	}
	printf ("\n");
}

int calc_design2 (int N, int L, int mode, int mode1, int mode2, int pnflag, int state[L][2 * N + 1], int i, int solution[N], int j)
{
	int k, l, NS, sum, flag;

	if (j == N)
	{
		for (k = 0; k < L; k++)
		{
			if (mode == 2 && k == 0)
				flag = 0;
			else if (mode2 == 0)
				flag = 0;
			else if (mode2 == 1)
			{
				if (state[k][N + i] == state[k][i])
					flag = 0;
				else
					flag = 1;
			}
			else if (mode2 == 2)
			{
				if (state[k][N + i] != state[k][2 * N])
					flag = 0;
				else
					flag = 1;
			}
			else if (mode2 == 3)
			{
				if (state[k][N + i] == state[k][i] && state[k][N + i] != state[k][2 * N])
					flag = 0;
				else
					flag = 1;
			}

			sum = 0;

			if (flag == 0)
			{
				for (l = 0; l < N; l++)
				{
					if (l == i)
						sum += solution[l];
					else
						sum += state[k][l] * solution[l];
				}
			}
			else
			{
				for (l = 0; l < N; l++)
				{
					if (solution[l] == 0)
						continue;
					else if (i == l)
						sum += solution[l];
					else if ((solution[l] > 0 && state[k][N + i] == 1) || (solution[l] < 0 && state[k][N + i] == 0))
					{
						sum += state[k][l] * state[k][N + l] * solution[l];
					}
					else if (state[k][l] + state[k][N + l] != 0)
					{
						sum += solution[l];
					}
				}
			}

			if ((sum > 0 && state[k][N + i] == 1) || (sum < 0 && state[k][N + i] == 0) || (sum == 0 && state[k][N + i] == state[k][i])) 
				continue;
			else
				return 0;
		}

		calc_links (N, solution);
		if (pnflag)
			print_solution (N, solution, i);

		return 1;
	}

	NS = 0;
	if (i != j)
	{
		solution[j] = -2;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		solution[j] = 0;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		solution[j] = 2;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
	}
	else
	{
		solution[j] = -1;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		solution[j] = 0;
		NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		if (mode1 == 1)
		{
			solution[j] = 1;
			NS += calc_design2 (N, L, mode, mode1, mode2, pnflag, state, i, solution, j + 1);
		}
	}

	return NS;
}

int calc_design (int N, int L, int mode, int mode1, int mode2, int pnflag, int pathway[L][N], int i)
{
	int j, k, state[L - 1][2 * N + 1], solution[N];

	for (j = 0; j < L - 1; j++)
	{
		for (k = 0; k < N; k++)
		{
			state[j][k] = pathway[j][k];
			state[j][N + k] = pathway[j + 1][k];
		}

		if (j != L - 2)
			state[j][2 * N] = pathway[j + 2][i];
		else
		{
			if (mode == 1)
				state[j][2 * N] = pathway[1][i];
			else
				state[j][2 * N] = pathway[j + 1][i];
		}
	}

	if (mode == 2)
	{
		if (i == 0)
			return 1;
		else
		{
			if (pathway[0][i] == pathway[1][i])
				solution[0] = 0;
			else if (pathway[1][i] == 1)
				solution[0] = 2;
			else if (pathway[1][i] == 0)
				solution[0] = -2;

			return calc_design2 (N, L - 1, mode, mode1, mode2, pnflag, state, i, solution, 1);
		}
	}
	else
	{
		return calc_design2 (N, L - 1, mode, mode1, mode2, pnflag, state, i, solution, 0);
	}
}

int main (int argc, char *argv[])
{
	int i, j, N, L, mode, mode1, mode2, pnflag;
	FILE *fp;
	char buf[256];

	if (argc == 1)
	{
		print_help ();
		exit (0);
	}

	mode = atoi (argv[1]);
	mode1 = atoi (argv[2]);
	mode2 = atoi (argv[3]);
	pnflag = atoi (argv[4]);
	N = atoi (argv[5]);
	L = atoi (argv[6]);

	int pathway[L][N], NS2[N], NS2_min;
	double NS;

	fp = fopen (argv[7], "r");
	while (fgets (buf, 256, fp))
	{
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				fscanf (fp, "%d", &pathway[i][j]);
			}
			fgets (buf, 256, fp);
		}

		NS = 1.0; minlinks = 0;
		for (i = 0; i < N; i++)
		{
			minlink = N;
			NS2[i] = calc_design (N, L, mode, mode1, mode2, pnflag, pathway, i);
			NS *= NS2[i];
			if (NS2[i] == 0)
				break;
			if (i == 0 || NS2[i] < NS2_min)
				NS2_min = NS2[i];
			minlinks += minlink;
		}

		if (NS == 0)
			continue;

		printf ("c %3d %4d %lg\n", minlinks, NS2_min, NS);
		for (i = 0; i < L; i++)
		{
			for (j = 0; j < N; j++)
			{
				printf ("%d ", pathway[i][j]);
			}
			printf ("\n");
		}
	}
	fclose (fp);
}
